home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_12_06 / prince / hidexpf.c < prev    next >
Text File  |  1994-01-24  |  1KB  |  34 lines
  1. /* _Expf function */
  2. #define _INCLUDE_XOPEN_SOURCE
  3. #include <math.h>
  4. #include "xmath.h"
  5.  
  6. #define hugexp (double)HUGE_EXP
  7.  
  8. int _Expf(double *px)
  9. {
  10.   /* Compute e^(*px), x finite */
  11.  
  12.   double y, g, even;
  13.   int xexp, inc;
  14.   const static double round[] = {.5, -.5};
  15. #if 0
  16. /* Works OK but slow on many machines */
  17.   xexp = *px * M_LOG2E + round[*px < 0];
  18. #else
  19.   /* VAX requires -(((short *)px)[1]>>15) */
  20.   xexp = *px * M_LOG2E + round[((unsigned int *) px)[_D0] >> 31];
  21. #endif
  22.   g = *px - xexp * M_LN2;
  23.   g *= (y = g * g) + 60.09114349;
  24.   even = y * 12.01517514 + 120.18228722;
  25.   *px = (even + g) / (even - g);
  26.   /* Limit exponent to wide enough range to cause over/underflow upon
  27.    * conversion to float without getting into double over/underflow */
  28.   inc = g > 0;
  29.   if (xexp > DBL_MAX_EXP - inc) xexp = DBL_MAX_EXP - inc;
  30.   if (xexp < DBL_MIN_EXP - inc) xexp = DBL_MIN_EXP - inc;
  31.   ((unsigned long *) px)[_D0] += xexp << _DOFF;    /* scale by 2^xexp */
  32.   return xexp + inc;        /* exponent = ceil(log2(*px)) */
  33. }
  34.